function [ModBFmat,QBFmat,modfn,phase_mod] = Beamforming(qbit,L,TRX,fc,BW_des,Ptx,showfig)
% qbit=1;

% Nsub=256;   % no. of subcarriers for rate computation
Ngrid = 1000;
% L = 15e-1;
% 
% factor = 0.3;
% x1 = 40e-1*factor;
% y1 = 0e-1*factor;
% z1 = 90e-1*factor;
% x2 = -18e-1*factor;
% y2 = 0e-1*factor;
% z2 = 20e-1*factor;
x1 = TRX(1);y1 = TRX(2);z1 = TRX(3);
x2 = TRX(4);y2 = TRX(5);z2 = TRX(6);
z12 = z1^2;z22 = z2^2;
% 
% fc = 30e9;
c = 3e8;
% BW_des = 1.5e9;
% Ptx = 0;%dBm


Ptx_lin = 10.^((Ptx-30)/10); % in watts
Lb = c/fc;
spac = Lb/2;
hp = 6.625e-34;
kb = 1.38e-23;
Tk = 290;


LIS_f = 0;
N = round(L/spac);

xcoord=linspace(-L/2,L/2,N);
ycoord=xcoord;
[xx,yy]=meshgrid(xcoord,ycoord);
coord_mat=xx+1i*yy;


alpha = atan((y1-y2)/(x1-x2));
xc = (x1+x2)/2;
yc = (y1+y2)/2;
xxp = (xx-xc)*cos(alpha)+(yy-yc)*sin(alpha);
yyp = -(xx-xc)*sin(alpha)+(yy-yc)*cos(alpha);



% L: x = xi + (xj-xi)t;
%    y = yi + (yj-yi)t;
Lxi = [xxp(1,1),xxp(1,end),xxp(end,end),xxp(end,1)];
Lyi = [yyp(1,1),yyp(1,end),yyp(end,end),yyp(end,1)];
Ldx = circshift(Lxi,-1)-Lxi;
Ldy = circshift(Lyi,-1)-Lyi;
L = [Lxi;Ldx;Lyi;Ldy];

u = (x1-xc)*cos(alpha)+(y1-yc)*sin(alpha);% x coordation of tx(x1,y1,z1)

a = (sqrt((xx-x1).^2+(yy-y1).^2+z1^2)+sqrt((xx-x2).^2+(yy-y2).^2+z2^2))/2;

a_vec = linspace(min(a(:)),max(a(:)),Ngrid)';% column vector
da_vec = a_vec(2)-a_vec(1);

b_vec = sqrt(a_vec.^2-u^2);
x0_vec = 1/4*u*(z22-z12)./b_vec.^2;
k_vec = sqrt(1+(x0_vec/u).^2-1/2*(z12+z22)./b_vec.^2);
A = b_vec.^2.*Ldx.^2+a_vec.^2.*Ldy.^2;
B = b_vec.^2 .*Ldx .*(Lxi-x0_vec)+a_vec.^2 .* Ldy .*Lyi;
C = b_vec.^2 .* (Lxi-x0_vec).^2 + a_vec.^2 .*Lyi.^2 -k_vec.^2.*a_vec.^2.*b_vec.^2;

t1 = (-B-sqrt(B.^2-A.*C))./A;
t2 = (-B+sqrt(B.^2-A.*C))./A;

isReal = (imag(t1)==0);
INrange1 = (t1<=1&t1>=0)&isReal;
INrange2 = (t2<=1&t2>=0)&isReal;
INrange = INrange1|INrange2;
cos_t1 = (Lxi-x0_vec+Ldx.*(t1))./k_vec./a_vec;
sin_t1 = (Lyi+Ldy.*t1)./k_vec./b_vec;
theta1 = mod(imag(log(cos_t1+1i*sin_t1)),2*pi);% bk sin(θ)= yi + (yt-yi)*t
cos_t2 = (Lxi-x0_vec+Ldx.*(t2))./k_vec./a_vec;
sin_t2 = (Lyi+Ldy.*t2)./k_vec./b_vec;
theta2 = mod(imag(log(cos_t2+1i*sin_t2)),2*pi);% bk sin(θ)= yi + (yt-yi)*t

dtheta1 = mod(theta1-circshift(theta2,1,2),2*pi).*INrange1.*circshift(INrange2,1,2);
dtheta2 = mod(theta1-circshift(theta2,2,2),2*pi).*INrange1.*circshift(INrange2,2,2).*circshift(~INrange,1,2);
dtheta3 = mod(theta1-circshift(theta2,3,2),2*pi).*INrange1.*circshift(INrange2,3,2).*circshift(~INrange,1,2).*circshift(~INrange,2,2);
dtheta4 = mod(theta1-circshift(theta2,4,2),2*pi).*INrange1.*circshift(INrange2,4,2).*circshift(~INrange,1,2).*circshift(~INrange,2,2).*circshift(~INrange,3,2);
t_proj = ((x0_vec-Lxi).*Ldx+(-Lyi).*Ldy)./(Ldx.^2+Ldy.^2);
inRIS = (sum(isReal,2)==0).*all(t_proj>0&t_proj<1,2);

modfn1 = sum(dtheta1,2)+sum(dtheta2,2)+sum(dtheta3,2)+sum(dtheta4,2)+inRIS*2*pi;
modfn1 = (a_vec.^2/2.*(-(z22-z12).^2/4./b_vec.^5+(z12+z22)./b_vec.^3)+k_vec.^2.*(b_vec.^2+1/2*u*u)./b_vec).*modfn1;

dstheta1 = (sin(theta1)-circshift(sin(theta2),1,2)).*INrange1.*circshift(INrange2,1,2);
dstheta2 = (sin(theta1)-circshift(sin(theta2),2,2)).*INrange1.*circshift(INrange2,2,2).*circshift(~INrange,1,2);
dstheta3 = (sin(theta1)-circshift(sin(theta2),3,2)).*INrange1.*circshift(INrange2,3,2).*circshift(~INrange,1,2).*circshift(~INrange,2,2);
dstheta4 = (sin(theta1)-circshift(sin(theta2),4,2)).*INrange1.*circshift(INrange2,4,2).*circshift(~INrange,1,2).*circshift(~INrange,2,2).*circshift(~INrange,3,2);

modfn2 = sum(dstheta1,2)+sum(dstheta2,2)+sum(dstheta3,2)+sum(dstheta4,2);
modfn2 = -a_vec*u.*(z22-z12)./(b_vec.^3).*k_vec.*modfn2;

dstheta1 = (sin(2*theta1)-circshift(sin(2*theta2),1,2)).*INrange1.*circshift(INrange2,1,2);
dstheta2 = (sin(2*theta1)-circshift(sin(2*theta2),2,2)).*INrange1.*circshift(INrange2,2,2).*circshift(~INrange,1,2);
dstheta3 = (sin(2*theta1)-circshift(sin(2*theta2),3,2)).*INrange1.*circshift(INrange2,3,2).*circshift(~INrange,1,2).*circshift(~INrange,2,2);
dstheta4 = (sin(2*theta1)-circshift(sin(2*theta2),4,2)).*INrange1.*circshift(INrange2,4,2).*circshift(~INrange,1,2).*circshift(~INrange,2,2).*circshift(~INrange,3,2);

modfn3 = sum(dstheta1,2)+sum(dstheta2,2)+sum(dstheta3,2)+sum(dstheta4,2);
modfn3 = -1/4.*k_vec.^2*u^2./b_vec.*modfn3;

modfn = modfn1+modfn2+modfn3;

if showfig
    close all
    [ccc,h] = contour(xcoord,ycoord,a,20);
    % clabel(ccc,h);
    figure;

    da = a_vec(end)-a_vec(1);
    a_vec_p = [linspace(a_vec(1)-da*0.5,a_vec(1),10)';a_vec;linspace(a_vec(end),a_vec(end)+da*0.5,10)'];
    modfn_p = [zeros(10,1);modfn;zeros(10,1)];    
    plot(a_vec_p*2,modfn_p)
    xlabel('RX到TX传输距离 l/m')
    ylabel('幅度|a(l)|')
    
    modfn_pappr = [zeros(10,1);mean(modfn).*ones(size(modfn));zeros(10,1)];
%     plot(a_vec_p*2/c*1e9,modfn_p,'LineWidth',1.5)
% xlabel('总传播时间\tau/ns')
% ylabel('调相波幅度|p(\tau)|')
end
Distmat = a*2;
distgrid = a_vec*2;
Distmat1 = sqrt((xx-x1).^2+(yy-y1).^2+z1^2);
Distmat2 = sqrt((xx-x2).^2+(yy-y2).^2+z2^2);

 %% RX parameters
%     Distmat_2D=abs(coord_mat-drx_x);
%     Distmat=sqrt(abs(Distmat_2D).^2+drx_z^2);%distance from tx points to rx
%     Dist_arr=abs(coord_mat);% distance between a point and the center.
    Hmat=(Lb/(2*pi))*exp(-1i*2*pi*fc*Distmat/c)./(Distmat1)./Distmat2;% channel response
%     Hmat(Dist_arr>Radius)=0;
    BFmat=exp(-1i*angle(Hmat));% traditional angle design
    QBFmat=quantz(exp(-1i*angle(Hmat)),qbit);% old distance matrix




P = cumsum(modfn.^2);
E = P(end);
phase_mod = 2*pi*BW_des*cumsum(P)/E*2*(a_vec(2)-a_vec(1))/c-pi*BW_des*2*(a_vec-a_vec(1))/c;


if showfig
modfn_2 = (2*modfn).^2;
freq_temp=cumsum(modfn_2)/Ngrid;% integration
scal=(2*pi*BW_des/c)/(freq_temp(end)-freq_temp(1));%  (39)
freq_temp=freq_temp*scal;  % (39)
freq_temp=freq_temp+(-pi*BW_des/c)-freq_temp(1); %(39)
phase_mod=cumsum(freq_temp)*(distgrid(2)-distgrid(1)); % as a function of rgrid- integration factor - (2*Radius)/Ngrid

plot(a_vec_p*2/c*1e9,modfn_p,'k','LineWidth',2)
hold on
grid on
plot(a_vec_p*2/c*1e9,modfn_pappr,'b--','LineWidth',2)
xlabel({'Propagation delay {\it t} (ns)'})
ylabel('Reflection intensity of Fresnel zone {\it v_t(t)}')
legend('Exact values','Approximate values','FontSize',16)
set(gca, 'FontSize', 16);
set(gcf, 'Position', [100 100 600 500]); % 调整当前图窗的大小
set(0, 'DefaultAxesFontName', 'Times New Roman');

subplot(2,2,1)
plot(a_vec_p*2/c*1e9,modfn_p,'LineWidth',1.5)
xlabel({'总传播时间\tau/ns','(a)'})
ylabel('调相波幅度|p(\tau)|')
subplot(2,2,2)
plot(a_vec*2/c*1e9, freq_temp*c,'LineWidth',1.5)
xlabel({'总传播时间\tau/ns','(b)'})
ylabel('瞬时角频率\Psi''(\tau)/rad\cdot s^{-1}')
subplot(2,2,3)

plot(a_vec*2/c*1e9, mod(phase_mod,2*pi),'LineWidth',1.5)
xlabel({'总传播时间\tau/ns','(c)'})
ylabel('设计相位\Psi(\tau)/rad')
subplot(2,2,4)
plot(a_vec*2/c*1e9,modfn.*cos(phase_mod),'LineWidth',1.5)
hold on
plot(a_vec*2/c*1e9,modfn.*sin(phase_mod),'--','LineWidth',1.5)
plot(90,14,'LineWidth',0.01)
hold off
xlabel({'总传播时间\tau/ns','(d)'})
ylabel('调频波波形')
legend('实部','虚部')
end
indx_valid=ones(size(a))>0;
GMod_off=exp(1i*compute_phase(Distmat,phase_mod,distgrid,indx_valid));% from ψ(u) to ψ(x,y)
ModBFmat=quantz(GMod_off.*BFmat,qbit);
ModBFmat=ModBFmat.*(indx_valid);
QBFmat=QBFmat.*(indx_valid);    %old 


        
        
        